1 Setup

1.1 Libraries used

Load libraries necessary to run the script.

library(tidyverse) # general
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0      ✔ purrr   1.0.1 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.5.0 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(ggdist) # plot distributions
library(ggside) # plot distributions on the plot margins
## Registered S3 method overwritten by 'ggside':
##   method from   
##   +.gg   ggplot2
library(easystats) # estimate and process statistical models
## # Attaching packages: easystats 0.6.0 (red = needs update)
## ✔ bayestestR  0.13.0   ✔ correlation 0.8.3 
## ✔ datawizard  0.6.5    ✖ effectsize  0.8.2 
## ✖ insight     0.18.8   ✖ modelbased  0.8.5 
## ✖ performance 0.10.1   ✖ parameters  0.20.1
## ✖ report      0.5.5    ✔ see         0.7.4 
## 
## Restart the R-Session and update packages in red with `easystats::easystats_update()`.
library(patchwork) # combine multiple plots
library(modelsummary) # for tables
## 
## Attaching package: 'modelsummary'
## 
## The following object is masked from 'package:parameters':
## 
##     supported_models
## 
## The following object is masked from 'package:insight':
## 
##     supported_models

1.2 Reusable functions

1.2.1 Plot distributions


plot_distributions <- function(df, col) {
  p1 <- df[c("Participant", col)] |>
    group_by(Participant) |>
    mean_qi(na.rm = TRUE) |>
    rename("x" = all_of(col)) |>
    mutate(Participant = fct_reorder(Participant, x)) |>
    ggplot(aes(x = x, y = Participant)) +
    ggdist::geom_pointinterval(aes(color = Participant, xmin = .lower, xmax = .upper), interval_size_range = c(0.03, 0.03), fatten_point = 1) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_discrete(expand = c(0, 0)) +
    scale_color_material_d(guide = "none") +
    theme_modern() +
    theme(axis.text.y = element_blank(), axis.title.x = element_blank()) +
    labs(y = "Participants") +
    ggside::geom_xsidedensity(aes(y = stat(density)), fill = "#9C27B0", color = NA) +
    ggside::geom_xsidedensity(data = df, aes_string(x = col, y = "stat(density)"), fill = "black", color = NA, alpha = 0.2) +
    ggside::theme_ggside_void()

  p2 <- df[c("Video", col)] |>
    group_by(Video) |>
    mean_qi(na.rm = TRUE) |>
    rename("x" = all_of(col)) |>
    mutate(Video = fct_reorder(Video, x)) |>
    ggplot(aes(x = x, y = Video)) +
    ggdist::geom_pointinterval(aes(color = Video, xmin = .lower, xmax = .upper), , interval_size_range = c(0.03, 0.03), fatten_point = 1) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_discrete(expand = c(0, 0)) +
    scale_color_viridis_d(guide = "none") +
    theme_modern() +
    theme(axis.text.y = element_blank(), axis.title.x = element_blank()) +
    labs(y = "Videos") +
    ggside::geom_xsidedensity(aes(y = stat(density)), fill = col_green, color = NA) +
    ggside::geom_xsidedensity(data = df, aes_string(x = col, y = "stat(density)"), fill = "black", color = NA, alpha = 0.2) +
    ggside::theme_ggside_void()

  wrap_elements(
    p1 + p2 +
      plot_annotation(
        title = col,
        theme = theme(plot.title = element_text(hjust = 0.5, face = "bold"))
      )
  )
}

2 Analysis

2.1 Load data

folder <- rstudioapi::selectDirectory(
  caption = "Select Directory",
  label = "Select"
)


df <- read.csv(paste0(folder, '/9_analysis/Online Study/data/data_pruned_for_analysis_2023-03-01.csv')) %>% 
  mutate(
    Typology = fct_relevel(video_primary_category, "urban", "green"),
    Valence = datawizard::rescale(Valence, to = c(1, 7), range = c(1, 9)),
    Arousal = datawizard::rescale(Arousal, to = c(1, 7), range = c(1, 9))
  ) %>% 
  rename("Participant" = Anon_ID )

length(unique(df$Participant))
## [1] 377

video_summary <- df %>% 
   group_by(Video) %>% 
  tally() %>% 
  arrange(n)  %>% 
  summarise(m = mean(n), s = sd(n), r = paste( range(n), collapse = " - "))

Each video was seen 23.16 times on average (SD = 1.876, range = 19 - 28).

2.2 Summary data per video (stimulus)

dfenv <- df %>% 
  group_by(video_name_df_keep, Typology) %>% 
  summarise_at(.vars = c( "Restoration", "WillingnessToWalk","Presence", "Width", "Structure", "Beauty","Interest", "Scenic", "Familiarity",  "Valence","Arousal", "Crowdedness" ), .funs = function(x) mean(x, na.rm = T)) %>% 
  ungroup() %>% 
  mutate(Typology = recode_factor(Typology, "urban" = "Urban", "green" = "Green" )) %>% 
  rename(`Willingness to walk` = WillingnessToWalk)

2.3 Participant

Here we only provide an overview to confirm we have loaded the correct amount of participants – as a sanity check. Participant characteristics are analysed and reported by a different file (Table 1).

dfsub <- df |>
  group_by(Participant) |>
  select(Participant, Age, Sex, Language, Crowd_Preference_Mean, starts_with("ipip"), SSA_Mean, SPS_Total_Mean, SIAS_Total_Mean, NSS, Upbringing_Environment, Current_Environment) |>
  slice(1)

nrow(dfsub) # should be 377
## [1] 377

A total of 377 participants (Mean age = 30.8, SD = 10.0, range: [18, 67]; Sex: 49.3% females, 50.7% males, 0.0% other) participants were included for analysis.

(dfsub |>
  ggplot(aes(x = Age)) +
  geom_density(fill = "orange") +
  geom_vline(xintercept = mean(dfsub$Age), color = "red", linetype = "dashed") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_modern()) |
  (dfsub |>
    ggplot(aes(x = Sex, fill = Sex)) +
    geom_bar(stat = "count") +
    scale_y_continuous(expand = c(0, 0)) +
    scale_fill_see_d(guide = "none") +
    theme_modern())

2.4 Descriptives of rerponses


library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows

df %>% 
  select(Restoration:Arousal) %>% 
  datasummary(All(.) ~  (`Unique (#)` = N) + (`Missing (%)` = PercentMissing) + Mean + SD, data =., 
              output = "latex_tabular", title = "Descriptive statistics of responses to the stimuli", 
              label = "response_descriptives")  %>% 
    save_kable(file = "output/response_descriptives.tex")
## Warning: To compile a LaTeX document with this table, the following commands must be placed in the document preamble:
## 
## \usepackage{booktabs}
## \usepackage{siunitx}
## \newcolumntype{d}{S[
##     input-open-uncertainty=,
##     input-close-uncertainty=,
##     parse-numbers = false,
##     table-align-text-pre=false,
##     table-align-text-post=false
##  ]}
## 
## To disable `siunitx` and prevent `modelsummary` from wrapping numeric entries in `\num{}`, call:
## 
## options("modelsummary_format_numeric_latex" = "plain")
##  This warning appears once per session.

3 Correlation plot

library(see)
library(correlation)

p_corrplot_u <- dfenv %>% 
  filter(Typology == "Urban") %>% 
  select(WLW = `Willingness to walk`, Crowdedness, Valence, Arousal, Restoration, Presence, Beauty, Interest, Structure, Width) %>%
  correlation(., bayesian = T) %>% 
  summary(redundant = FALSE) %>%
  plot() + scale_x_discrete(position = "top") + labs(title = "Urban spaces") +
  theme_modern()
## Warning in genhypergeo_series_pos(U = c((n - 1)/2, (n - 1)/2), L = ((n + :
## Series not converged.

## Warning in genhypergeo_series_pos(U = c((n - 1)/2, (n - 1)/2), L = ((n + :
## Series not converged.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

p_corrplot_g <- dfenv %>% 
  filter(Typology == "Green") %>% 
  select(WLW = `Willingness to walk`, Crowdedness, Valence, Arousal, Restoration, Presence, Beauty, Interest, Structure, Width) %>%
  correlation(., bayesian = T) %>% 
  summary(redundant = FALSE) %>% 
  plot() + scale_x_discrete(position = "top") + labs(title = "Green spaces") +
  theme_modern()
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

p_corrplot_g



p_corr <- p_corrplot_u / p_corrplot_g + plot_annotation(title = "Correlation Matrix")
ggsave(p_corr, filename = "output/plots_corr.jpg", width = 13, height = 7, dpi = 300)

3.1 Video Ratings

3.1.1 Variable Distribution

3.2 Determinants of Perceived Crowdedness

3.2.1 Models

rm(preds, params, mods, vars)
## Warning in rm(preds, params, mods, vars): object 'preds' not found
## Warning in rm(preds, params, mods, vars): object 'params' not found
## Warning in rm(preds, params, mods, vars): object 'mods' not found
## Warning in rm(preds, params, mods, vars): object 'vars' not found
preds <- data.frame()
params <- data.frame()
mods <- list()
outcome_var <- "Crowdedness"
vars <- c("Beauty", "Scenic","Width", "Structure") # replaced familiarity with Beauty
i=1
for (i in 1:length(vars)) {
  var <- vars[i]
  print(var)

  f <- paste0(outcome_var, " ~ log1p(video_mean_pedcounts) * ", var, " * Typology + Familiarity  + (", var, "|Participant) + (1|Video)")

  model <- glmmTMB::glmmTMB(as.formula(f), data = df)
  mods[[i]] <- model
  
  # Parameters
  param <- parameters::parameters(model, effects = "fixed")
  param$Variable <- var
  params <- rbind(params, param)

  # Predicted
  pred <- estimate_relation(model, at = c("video_mean_pedcounts", var, "Typology"), length = c(50, 7))
  names(pred)[names(pred) == var] <- "Value"
  pred$Variable <- var
  pred$Label <- paste0("Main effect: ", format_p(param$p[3]), "\nInteraction: ", format_p(param$p[4]))

  preds <- rbind(preds, pred)
}
## [1] "Beauty"
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.4.1
## Current Matrix version is 1.5.1
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
## [1] "Scenic"
## [1] "Width"
## [1] "Structure"

3.2.2 plots

rm(p_crowd)
## Warning in rm(p_crowd): object 'p_crowd' not found
p_crowd <- preds |>
  filter(Predicted < 7) |>
  mutate(Value = as.factor(Value)) |>
  ggplot(aes(x = video_mean_pedcounts, y = Predicted)) +
  stat_density_2d(data = df, aes(y = Crowdedness, fill = ..density..), geom = "raster", contour = FALSE) +
  geom_line(aes(color = Typology, size = Value, group = interaction(Typology, Value))) +
  scale_fill_gradientn(colours = c("white", "grey"), guide = "none") +
  scale_size_manual(values = seq(from = .05, to = 2, length.out= 7), breaks= c(1:7)) +
  # scale_colour_viridis_b(option = "A") +
  # scale_x_log10(expand = c(0, 0), breaks = seq(2, 20, by = 2)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), breaks = c(0, 1, 3, 5, 7)) +
  labs(x = "Number of Pedestrians", y = "Perceived Crowdedness") +
  # coord_cartesian(ylim = c(1, 7)) +
  labs(alpha = "Rating", title = "Predictors of Crowdedness") +
  theme_modern(axis.title.space = 10) +
  theme(
    strip.placement = "outside",
    strip.text.x = element_text(hjust = 0, size = 12),
    axis.title.y = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5)
  ) +
  ggnewscale::new_scale_fill() +
  scale_color_manual(values = c("urban" = col_urban, "green" = col_green )) + # col_green
  scale_fill_manual(values = c("urban" = col_urban, "green" = col_green)) +
  ggside::geom_xsidedensity(data = df, aes(y = after_stat(scaled), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
  ggside::geom_ysidedensity(data = df |>
    select(all_of(c("Typology", vars))) |>
    pivot_longer(-Typology, values_to = "Predicted", names_to = "Variable"), aes(x = after_stat(scaled), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
  ggside::theme_ggside_void() +
  facet_wrap(~Variable, strip.position = "top", ncol = 2)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

p_crowd
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## Warning: Removed 12 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).

3.2.3 tables

library(modelsummary)
options("modelsummary_format_numeric_latex" = "plain")

names(mods) <- vars
modelsummary(mods, escape = FALSE,
  stars = TRUE, statistic = "[{conf.low} - {conf.high}]",
  coef_omit = "SD|Cor", 
  coef_rename = c(
    "log1p(mean_pedcounts)" = "Ped. Counts",
    "Typologygreen" = "Green",
    
    "StructureTypologygreen" = "Structure x Green",
    "BeautyTypologygreen" = "Beauty x Green",
    "WidthTypologygreen" = "Width x Green",
    "ScenicTypologygreen" = "Scenic x Green",
    
    "log1p(mean_pedcounts)Typologygreen"= "Ped. Counts x Green",
    "log1p(mean_pedcounts)Beauty" = "Ped. Counts x Beauty",
    "log1p(mean_pedcounts)Scenic" = "Ped. Counts x Scenic",
    "log1p(mean_pedcounts)Width" = "Ped. Counts x Width",
    "log1p(mean_pedcounts)Structure" = "Ped. Counts x Structure",

    "log1p(mean_pedcounts)ScenicTypologygreen" = "Ped. Counts x Scenic x Green",
    "log1p(mean_pedcounts)WidthTypologygreen" = "Ped. Counts x Width x Green",
    "log1p(mean_pedcounts)StructureTypologygreen" = "Ped. Counts x Structure x Green",
    "log1p(mean_pedcounts)BeautyTypologygreen" = "Ped. Counts x Beauty x Green"
  ),
  title = "Mixed effects models of Perceived Crowdedness \\label{tab:models_perceived_crowds}", 
  notes = "Ped. Counts : Log-transformed, average pedestrian count per stimulus",
  output = "output/table_models_perceived_crowds.tex"
)

3.3 Determinants of Beauty

rm(preds, mods, vars)
preds <- data.frame()
dat <- data.frame()
mods <- list()
vars <- c("Crowdedness", "Structure", "Scenic", "Width") #  "Interest",

for (i in 1:length(vars)) {
  var <- vars[i]
  print(var)

  f <- paste0("Beauty ~ poly(", var, ", 2) * Typology + Familiarity + (", var, "|Participant)")

  model <- glmmTMB::glmmTMB(as.formula(f), data = na.omit(df[c(var, "Beauty", "Typology", "Familiarity", "Participant")]))
  
  mods[[i]] <- model
  # Individual-level correlations
  # dfsub <- modelbased::estimate_grouplevel(model) |>
  #   modelbased::reshape_grouplevel() |>
  #   select(c(1, 3)) |>
  #   group_by(Participant) |>
  #   slice(1) |>
  #   data_rename(pattern = paste0("Participant_Coefficient_", var), paste0("Beauty_", var), regex = TRUE) |>
  #   merge(dfsub, by = "Participant")

  # Parameters
  param <- parameters::parameters(model, effects = "fixed")

  # Predicted
  pred <- estimate_relation(model, at = c(var, "Typology"), length = c(30))
  names(pred)[names(pred) == var] <- "Value"
  pred$Variable <- var
  pred$Label <- paste0(
    "Objective Crowding: ", format_p(param$p[2]),
    "\n", var, ": ", format_p(param$p[3]),
    "\nInteraction: ", format_p(param$p[4])
  )

  # Data for density
  dat <- rbind(dat, data.frame(Value = df[[var]], Predicted = df$Beauty, Variable = var))
  preds <- rbind(preds, pred)
}
## [1] "Crowdedness"
## [1] "Structure"
## [1] "Scenic"
## [1] "Width"

3.3.1 plots

p_beauty <- preds |>
  ggplot(aes(x = Value, y = Predicted)) +
  stat_density_2d(data = dat, aes(fill = ..density..), geom = "raster", contour = FALSE) +
  geom_line(aes(color = Typology), size = 1) +
  scale_fill_gradientn(colours = c("white", "grey"), guide = "none") +
  scale_color_manual(values = c("urban" = col_urban, "green" = col_green)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(y = "Beauty", title = "Predictors of Beauty") +
  theme_modern(axis.title.space = 10) +
  theme(
    axis.title.x = element_blank(),
    strip.placement = "outside",
    strip.text.x = element_text(hjust = 0.5, size = 12, face = "plain"),
    axis.title.y = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5)
  ) +
  ggnewscale::new_scale_fill() +
  scale_fill_manual(values = c("urban" = col_urban, "green" = col_green)) +
  ggside::geom_ysidedensity(data = mutate(df, Predicted = Beauty), aes(x = after_stat(scaled), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
  ggside::geom_xsidedensity(data = df |>
    select(all_of(c("Typology", vars))) |>
    pivot_longer(-Typology, values_to = "Value", names_to = "Variable"), aes(y = after_stat(scaled), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
  ggside::theme_ggside_void() +
  facet_wrap(~Variable, strip.position = "top", scales = "free_x")

p_beauty
## Warning: Removed 34 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 12 rows containing non-finite values (`stat_density()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).

3.3.2 Tables

names(mods) <- vars
modelsummary(mods, stars = TRUE)
Crowdedness Structure Scenic Width
(Intercept) 3.861*** 4.006*** 4.469*** 4.071***
(0.049) (0.045) (0.024) (0.045)
poly(Crowdedness, 2)1 −0.972
(2.226)
poly(Crowdedness, 2)2 −17.489***
(1.568)
Typologygreen 1.328*** 1.087*** 0.243*** 1.178***
(0.042) (0.040) (0.028) (0.036)
Familiarity 0.133*** 0.113*** 0.033*** 0.099***
(0.008) (0.008) (0.005) (0.008)
poly(Crowdedness, 2)1 × Typologygreen −18.689***
(4.482)
poly(Crowdedness, 2)2 × Typologygreen 9.603*
(4.064)
poly(Structure, 2)1 55.409***
(2.156)
poly(Structure, 2)2 4.805**
(1.542)
poly(Structure, 2)1 × Typologygreen 0.514
(4.849)
poly(Structure, 2)2 × Typologygreen −3.159
(4.236)
poly(Scenic, 2)1 137.104***
(1.270)
poly(Scenic, 2)2 −2.104*
(0.963)
poly(Scenic, 2)1 × Typologygreen 1.391
(3.208)
poly(Scenic, 2)2 × Typologygreen −4.448+
(2.626)
poly(Width, 2)1 59.713***
(2.295)
poly(Width, 2)2 5.092***
(1.534)
poly(Width, 2)1 × Typologygreen −17.449***
(3.576)
poly(Width, 2)2 × Typologygreen 1.021
(3.480)
SD (Intercept Participant) 0.777 0.998 0.536 1.041
SD (Crowdedness Participant) 0.155
Cor (Intercept~Crowdedness Participant) −0.614
SD (Structure Participant) 0.166
Cor (Intercept~Structure Participant) −0.848
SD (Scenic Participant) 0.091
Cor (Intercept~Scenic Participant) −0.917
SD (Width Participant) 0.178
Cor (Intercept~Width Participant) −0.857
SD (Observations) 1.290 1.216 0.728 1.201
Num.Obs. 11107 11102 11104 11105
R2 Marg. 0.136 0.195 0.688 0.193
R2 Cond. 0.366 0.519 0.797 0.539
AIC 38293.3 36891.5 25248.7 36703.2
BIC 38373.8 36972.0 25329.1 36783.6
ICC 0.3 0.4 0.4 0.4
RMSE 1.26 1.19 0.71 1.17
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
modelsummary(mods, escape = FALSE,
  stars = TRUE, statistic = "std.error", #  "[{conf.low} - {conf.high}]",
  coef_omit = "SD|Cor",
  coef_rename = c(
                  "poly(Crowdedness, 2)1" = "Crowdedness",
                  "poly(Crowdedness, 2)2" = "Crowdedness<sup> 2</sup>",
                  "Typologygreen" = "Green",
                  "Familiarity"  = "Familiarity",
                  
                  "poly(Crowdedness, 2)1Typologygreen" = "Crowdedness x Green",
                  "poly(Crowdedness, 2)2Typologygreen" = "Crowdedness<sup> 2</sup> x Green",
                  
                  "poly(Structure, 2)1" = "Structure",
                  "poly(Structure, 2)2" = "Structure<sup> 2</sup>",
                  
                  "poly(Structure, 2)1Typologygreen" = "Structure x Green",
                  "poly(Structure, 2)2Typologygreen" = "Structure<sup> 2</sup> x Green",
                  
                  "poly(Scenic, 2)1" = "Scenic",
                  "poly(Scenic, 2)2" = "Scenic<sup> 2</sup>",
                  
                  "poly(Scenic, 2)1Typologygreen" = "Scenic x Green",
                  "poly(Scenic, 2)2Typologygreen" = "Scenic<sup> 2</sup> x Green",
                  
                  "poly(Width, 2)1Typologygreen" = "Width x Green",
                  "poly(Width, 2)2Typologygreen" = "Width<sup> 2</sup> x Green"
                  
                  ),
  title = "Non-linear mixed effects models of Perceived Beauty \\label{tab:models_beauty}",
  # notes = "PC (log) : Average pedestrian count per stimulus" ,
  output = "output/table_models_perceived_beauty.tex"
)

3.4 Determinants of Psychological Reaction

rm(preds, dat, mods, vars)
preds <- data.frame()
dat <- data.frame()


mods <- list()
outcome_vars <- c("WillingnessToWalk", "Valence", "Arousal", "Restoration",  "Presence")

for (i in 1:length(outcome_vars)) {
  resp <- outcome_vars[i]
  print(resp)

  f <- paste(resp, " ~ Typology / (poly(Crowdedness, 2) * Beauty) + Familiarity + (1|Participant) + (1|Video)")
  f
  model <- glmmTMB::glmmTMB(as.formula(f), data = na.omit(df[c(resp, "Beauty", "Typology", "Crowdedness", "Familiarity", "Participant", "Video")]))
  mods[[i]] <- model
  param <- parameters::parameters(model, effects = "fixed")

  # Predictions
  preds <- estimate_relation(model, at = c("Typology", "Crowdedness", "Beauty"), length = c(50, 7)) |>
    mutate(Outcome = resp) |>
    rbind(preds)

  # Data for density
  dat <- rbind(
    dat,
    data.frame(
      Typology = df$Typology,
      Crowdedness = df$Crowdedness,
      video_mean_pedcounts = df$video_mean_pedcounts,
      Predicted = df[[resp]],
      Outcome = resp
    )
  )

  # Random (simplified model)
  # f <- paste(resp, " ~ mean_pedcounts + Beauty + (mean_pedcounts + Beauty |Participant)")
  # model <- glmmTMB::glmmTMB(as.formula(f), data = df)
  # # dfsub <- modelbased::estimate_grouplevel(model) |>
  #   modelbased::reshape_grouplevel(indices = "Coefficient") |>
  #   select(c(1, 3, 4)) |>
  #   group_by(Participant) |>
  #   slice(1) |>
  #   data_addprefix(paste0(resp, "_"), exclude = "Participant") |>
  #   merge(dfsub, by = "Participant")
}
## [1] "WillingnessToWalk"
## [1] "Valence"
## [1] "Arousal"
## [1] "Restoration"
## [1] "Presence"



names(mods) <- c("Willingness to Walk", "Valence", "Arousal", "Restoration",  "Presence")

3.4.1 plots psych


dat$Outcome <- recode(dat$Outcome, "WillingnessToWalk"="Willing.-to-walk")

p_psych <- preds |>
  dplyr::filter(Predicted <= 7) |>
  mutate(Beauty = fct_rev(as.factor(Beauty)),
         Outcome = recode_factor(Outcome, "WillingnessToWalk"="Willing.-to-walk"),
         Outcome = fct_relevel(Outcome, c("Willing.-to-walk", "Valence", "Arousal", "Restoration", "Interest", "Presence"))) |>
  ggplot(aes(x = Crowdedness, y = Predicted)) +
  stat_density_2d(data = dat, aes(fill = ..density..), geom = "raster", contour = FALSE) +
  geom_line(aes(size = Beauty, color = Typology, group = interaction(Beauty, Typology))) +
  scale_size_manual(values = seq(from = .05, to = 2, length.out= 7), breaks= c(1:7)) +
  scale_fill_gradientn(colours = c("white", "grey"), guide = "none") +
  scale_color_manual(values = c("urban" = col_urban, "green" = col_green)) +
  scale_alpha_ordinal(range = c(1, 0.1)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), breaks = c(1:7)) +
  theme_modern(axis.title.space = 10) +
  theme(
    strip.placement = "outside",
    strip.text.x = element_text(hjust = 0, size = 12),
    axis.title.y = element_text(face = "bold"),
    plot.title = element_text(face = "bold", hjust = 0.5)
  ) +
  labs(x = "Subjective Crowdedness", y = "Psychological Outcome", title = "Psychological effect of Crowdedness and its modulation by Beauty") +
  ggnewscale::new_scale_fill() +
  scale_fill_manual(values = c("urban" = col_urban, "green" = col_green)) +
  ggside::geom_xsidedensity(data=dat, aes(y=stat(density), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
  ggside::geom_ysidedensity(data=dat, aes(x=stat(density), fill = Typology), color = NA, alpha = 0.5, adjust = 2) +
  ggside::theme_ggside_void() +
  facet_wrap(~Outcome, strip.position = "left", scales = "free")

p_psych

3.4.2 Table output

modelsummary(mods, stars = TRUE)
Willingness to Walk  Valence Arousal Restoration Presence
(Intercept) 0.078 1.675*** 2.893*** −0.320*** 2.490***
(0.053) (0.044) (0.067) (0.064) (0.070)
Typologygreen 0.064 0.028 0.007 0.540** 0.062
(0.155) (0.127) (0.196) (0.173) (0.138)
Familiarity 0.052*** 0.066*** 0.013 0.013+ 0.151***
(0.007) (0.006) (0.009) (0.008) (0.006)
Typologyurban × poly(Crowdedness, 2)1 −5.507 −8.341** 12.964** 0.898 15.066***
(3.460) (2.839) (4.398) (3.908) (3.062)
Typologygreen × poly(Crowdedness, 2)1 −17.504 −36.489*** −1.308 −19.692 15.132
(13.306) (10.936) (16.952) (15.026) (12.081)
Typologyurban × poly(Crowdedness, 2)2 −26.159*** −29.599*** 12.199** −16.120*** 3.829
(3.361) (2.753) (4.265) (3.777) (3.071)
Typologygreen × poly(Crowdedness, 2)2 −25.949* −40.108*** −3.835 −22.872+ −16.571
(12.189) (9.979) (15.459) (13.692) (11.337)
Typologyurban × Beauty 0.758*** 0.638*** 0.206*** 0.799*** 0.255***
(0.008) (0.007) (0.010) (0.010) (0.007)
Typologygreen × Beauty 0.741*** 0.636*** 0.083* 0.763*** 0.250***
(0.026) (0.021) (0.032) (0.028) (0.023)
Typologyurban × poly(Crowdedness, 2)1 × Beauty 0.625 0.952 6.878*** −1.104 −2.505***
(0.731) (0.599) (0.930) (0.824) (0.662)
Typologygreen × poly(Crowdedness, 2)1 × Beauty 0.292 4.187* 4.824+ 1.049 −2.959
(2.267) (1.861) (2.883) (2.552) (2.087)
Typologyurban × poly(Crowdedness, 2)2 × Beauty 3.214*** 4.649*** −3.547*** 1.357+ −0.503
(0.723) (0.591) (0.917) (0.811) (0.667)
Typologygreen × poly(Crowdedness, 2)2 × Beauty 2.136 5.612*** −0.186 2.660 1.931
(2.072) (1.696) (2.629) (2.327) (1.931)
SD (Intercept Participant) 0.528 0.461 0.670 0.742 1.107
SD (Intercept Video) 0.175 0.151 0.229 0.215 0.000
SD (Observations) 1.010 0.825 1.280 1.129 0.951
Num.Obs. 10899 11009 11009 11107 11106
R2 Marg. 0.537 0.554 0.119 0.486 0.239
R2 Cond. 0.645 0.669 0.326 0.650
AIC 32248.3 28168.2 37791.7 35523.6 31831.3
BIC 32365.0 28285.1 37908.6 35640.7 31948.3
ICC 0.2 0.3 0.2 0.3
RMSE 0.99 0.80 1.25 1.10 0.94
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

modelsummary(mods,
  stars = TRUE, statistic = "[{conf.low} - {conf.high}]", escape = T,
  coef_omit = "SD|Cor", 
  coef_rename = c("Typologygreen"="Green",
                  "Familiarity"  = "Familiarity",
                  "Typologyurbanpoly(Crowdedness, 2)1"="Urban * Crowdedness",
                  "Typologygreenpoly(Crowdedness, 2)1"="Green * Crowdedness",
                  "Typologyurbanpoly(Crowdedness, 2)2"="Urban * Crowdedness<sup> 2</sup>",
                  "Typologygreenpoly(Crowdedness, 2)2"="Green * Crowdedness<sup> 2</sup>",
                  "TypologyurbanBeauty" = "Urban * Beauty",
                  "TypologygreenBeauty" = "Green * Beauty",
                  "Typologyurbanpoly(Crowdedness, 2)1Beauty" = "Urban * Crowdedness * Beauty",
                  "Typologygreenpoly(Crowdedness, 2)1Beauty" = "Green * Crowdedness * Beauty",
                  "Typologyurbanpoly(Crowdedness, 2)2Beauty" = "Urban * Crowdedness<sup> 2</sup> * Beauty",
                  "Typologygreenpoly(Crowdedness, 2)2Beauty" = "Green * Crowdedness<sup> 2</sup> * Beauty"
                  ),
  title = "Non-linear mixed effects models of psychological reaction\\label{tab:models_reaction}",
  notes = "PC (log) : Average pedestrian count per stimulus" ,
  output = "output/table_models_psych_reaction.tex"
)

3.5 joint plots


p_reg <- ((p_crowd + guides(fill = "none", color = "none") | p_beauty + theme(legend.position = "none")) / p_psych + theme(legend.position = 'bottom') )  +plot_annotation(tag_levels = "A")

p_reg
## Warning: Removed 12 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
## Warning: Removed 34 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 12 rows containing non-finite values (`stat_density()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
## Warning: Removed 420 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 15 rows containing non-finite values (`stat_density()`).
## Warning: Removed 405 rows containing non-finite values (`stat_density()`).
## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.


ggsave(p_reg, filename = "output/plots_reg.jpg", width = 13, height = 13, dpi = 300)
## Warning: Removed 12 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
## Warning: Removed 34 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 12 rows containing non-finite values (`stat_density()`).
## Warning: Removed 22 rows containing non-finite values (`stat_density()`).
## Warning: Removed 420 rows containing non-finite values (`new_stat_density2d()`).
## Warning: Removed 15 rows containing non-finite values (`stat_density()`).
## Warning: Removed 405 rows containing non-finite values (`stat_density()`).
## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
pvalue <- function(x, ...) {
    # Construct vectors of data y, and groups (strata) g
    y <- unlist(x)
    g <- factor(rep(1:length(x), times=sapply(x, length)))
    if (is.numeric(y)) {
        # For numeric variables, perform a standard 2-sample t-test
        # p <- t.test(y ~ g)$p.value
        p <- BayesFactor::ttestBF(x=y, y= g)
        p <- p@bayesFactor$bf
            # print(p)
    } else {
        # For categorical variables, perform a chi-squared test of independence
        p <- chisq.test(table(y, g))$p.value
    }
    # Format the p-value, using an HTML entity for the less-than sign.
    # The initial empty string places the output on the line below the variable label.
    c("", sub("<", "&lt;", format.pval(p, digits=3, eps=0.001)))
}

library(kableExtra)

t1 <- table1::table1(~ Valence + Arousal + Restoration + Beauty + `Willingness to walk` + Interest + Familiarity + Structure + Presence  + Crowdedness  | Typology, overall = F,
    data=dfenv, extra.col=list(`BF`=pvalue))
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
## t is large; approximation invoked.
t1
Urban
(N=408)
Green
(N=72)
BF
Valence
Mean (SD) 4.75 (0.638) 5.81 (0.525) 1175
Median [Min, Max] 4.81 [2.53, 6.19] 5.94 [3.36, 6.71]
Arousal
Mean (SD) 3.90 (0.547) 3.29 (0.546) 1018
Median [Min, Max] 3.91 [2.14, 5.33] 3.33 [1.97, 4.46]
Restoration
Mean (SD) 3.25 (0.817) 4.85 (0.532) 606
Median [Min, Max] 3.23 [0.870, 5.36] 4.95 [2.91, 5.74]
Beauty
Mean (SD) 4.41 (0.761) 5.89 (0.514) 958
Median [Min, Max] 4.36 [1.96, 6.41] 6.02 [4.27, 6.83]
Willingness to walk
Mean (SD) 3.63 (0.750) 4.82 (0.584) 787
Median [Min, Max] 3.68 [1.30, 5.44] 4.98 [2.52, 5.74]
Interest
Mean (SD) 4.33 (0.832) 5.35 (0.608) 939
Median [Min, Max] 4.33 [2.22, 6.24] 5.40 [3.80, 6.23]
Familiarity
Mean (SD) 4.10 (0.559) 4.49 (0.602) 1130
Median [Min, Max] 4.13 [2.05, 5.68] 4.52 [2.82, 5.76]
Structure
Mean (SD) 4.37 (0.989) 5.48 (0.606) 847
Median [Min, Max] 4.45 [1.40, 6.54] 5.61 [3.29, 6.32]
Presence
Mean (SD) 4.23 (0.332) 4.71 (0.339) 1419
Median [Min, Max] 4.25 [2.91, 5.13] 4.77 [3.82, 5.48]
Crowdedness
Mean (SD) 4.42 (1.42) 2.99 (1.47) 509
Median [Min, Max] 4.56 [1.22, 6.88] 2.84 [1.14, 5.96]
t1 %>% data.frame() %>% filter(str_detect(X. , "Median", negate = T)) %>% rename(" " = "X.") %>% kbl(., format = "latex", escape = T, booktabs = T, label = "stimuli_descriptives", caption = "Descriptive statistics of emotional reaction and environmental appraisals for stimuli, split by typology") %>% row_spec(seq(2, 20, 2), bold=T) %>% add_indent(seq(3, 21, 2), ) %>% footnote(number = "Bayesian t-tests; reporting the bayes factor.") %>% save_kable(file = "output/stimuli_table1.tex", )
library(ggdist)

penv <- dfenv %>% 
  pivot_longer(Restoration : Crowdedness) %>% 
  ggplot(aes(x = name, y = value, colour = Typology)) +
  # geom_boxplot( position = position_dodge()) +
  # geom_line(aes(group = name_df_keep),  alpha =.2)+
  geom_jitter(position = position_jitterdodge(), alpha = 0.5) +
  ggdist::stat_halfeye(position = position_dodge(), aes(fill = Typology),
    ## custom bandwidth
    adjust = .5,
    ## adjust height
    width = .6,
    ## move geom to the right
    justification = -.4,
    ## remove slab interval
    .width = 0,
  point_colour = NA) +
  scale_color_manual(values = c("Urban" = col_urban, "Green" = col_green)) +
  scale_fill_manual(values = c("Urban" = col_urban, "Green" = col_green)) +
  theme_classic() +
  # gghalves
  labs(y = "Rating", x = "Dimension")

penv

ggsave(penv, filename = "output/plots_envstats.jpg", width = 13, height = 13/2, dpi = 300)
library(see) # for plotting
library(ggraph) # needs to be loaded

names(dfenv) <- make.names(names(dfenv))
cor_urb <- dfenv %>% 
  filter(Typology == "Urban") %>% 
  ungroup() %>% 
  select(Willingness.to.walk, Restoration, Interest, Valence, Arousal, Crowdedness, Scenic, Width ) %>% 
  correlation(partial = TRUE, standardize_names = F) %>% 
  filter(p < 0.01) %>%
  plot()
cor_urb


cor_green <- dfenv %>% 
  filter(Typology == "Green") %>% 
  ungroup() %>% 
  select(Willingness.to.walk, Restoration, Interest, Valence, Arousal, Crowdedness, Scenic, Width ) %>% 
  correlation(partial = TRUE) %>% 
  # filter(p < 0.01) %>% 
  plot()
cor_green

3.6 Pics


dfenv %>% 
  ggplot(aes(Crowdedness,Beauty,  colour = Typology)) +
  geom_jitter(alpha = 0.5, position = position_jitterdodge(jitter.height = .3, jitter.width = .3))+ 
  geom_smooth(method = "gam" ) 
## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'
## Warning: `position_jitterdodge()` requires non-overlapping x intervals

  scale_color_manual(values = c("urban" = col_urban, "green" = col_green)) 
## <ggproto object: Class ScaleDiscrete, Scale, gg>
##     aesthetics: colour
##     axis_order: function
##     break_info: function
##     break_positions: function
##     breaks: waiver
##     call: call
##     clone: function
##     dimension: function
##     drop: TRUE
##     expand: waiver
##     get_breaks: function
##     get_breaks_minor: function
##     get_labels: function
##     get_limits: function
##     guide: legend
##     is_discrete: function
##     is_empty: function
##     labels: waiver
##     limits: function
##     make_sec_title: function
##     make_title: function
##     map: function
##     map_df: function
##     n.breaks.cache: NULL
##     na.translate: TRUE
##     na.value: grey50
##     name: waiver
##     palette: function
##     palette.cache: NULL
##     position: left
##     range: <ggproto object: Class RangeDiscrete, Range, gg>
##         range: NULL
##         reset: function
##         train: function
##         super:  <ggproto object: Class RangeDiscrete, Range, gg>
##     rescale: function
##     reset: function
##     scale_name: manual
##     train: function
##     train_df: function
##     transform: function
##     transform_df: function
##     super:  <ggproto object: Class ScaleDiscrete, Scale, gg>

dfenv %>% 
  group_by(Typology) %>% 
  arrange(Beauty) %>% 
  slice(1:5, (n()-5): n()) %>% 
  select(video_name_df_keep, Beauty)
## Adding missing grouping variables: `Typology`
## # A tibble: 22 × 3
## # Groups:   Typology [2]
##    Typology video_name_df_keep                                            Beauty
##    <fct>    <chr>                                                          <dbl>
##  1 Urban    Walking in USTI NAD LABEM _ Czech Republic - 4K 60fps (UHD).…   1.96
##  2 Urban    Walking in LIVERPOOL _ UK - 4K 60fps (UHD).mp4_17.mp4           2.45
##  3 Urban    4K Virtual walking Vienna - Augarten Park (1).mp4_3.mp4         2.64
##  4 Urban    Walking in GEORGETOWN _ Penang (Malaysia)  - Jetties and Old…   2.81
##  5 Urban    Walking in GEORGETOWN _ Penang (Malaysia)  - Jetties and Old…   2.87
##  6 Urban    4K Virtual walking Vienna - Augarten Park (1).mp4_18.mp4        6.10
##  7 Urban    Walking in SINGAPORE - Chinatown to Marina Bay and Gardens -…   6.16
##  8 Urban    Walking in SALZBURG _ Austria - Bells in the Old Town (2019)…   6.19
##  9 Urban    Walking in LIVERPOOL _ UK - 4K 60fps (UHD).mp4_10.mp4           6.25
## 10 Urban    Walking in ZURICH _ Switzerland - 4K 60fps (UHD).mp4_2.mp4      6.28
## # … with 12 more rows

3.7 Session information

sessionInfo()
## R version 4.2.2 (2022-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.0
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggraph_2.1.0       kableExtra_1.3.4   modelsummary_1.3.0 patchwork_1.1.2   
##  [5] see_0.7.4          report_0.5.5       parameters_0.20.1  performance_0.10.1
##  [9] modelbased_0.8.5   insight_0.18.8     effectsize_0.8.2   datawizard_0.6.5  
## [13] correlation_0.8.3  bayestestR_0.13.0  easystats_0.6.0    ggside_0.2.2      
## [17] ggdist_3.2.1       forcats_0.5.2      stringr_1.5.0      dplyr_1.0.10      
## [21] purrr_1.0.1        readr_2.1.3        tidyr_1.2.1        tibble_3.1.8      
## [25] ggplot2_3.4.0      tidyverse_1.3.2   
## 
## loaded via a namespace (and not attached):
##   [1] googledrive_2.0.0      ggnewscale_0.4.8       minqa_1.2.5           
##   [4] colorspace_2.0-3       ellipsis_0.3.2         estimability_1.4.1    
##   [7] fs_1.5.2               rstudioapi_0.14        glmmTMB_1.1.5         
##  [10] farver_2.1.1           graphlayouts_0.8.4     MatrixModels_0.5-1    
##  [13] ggrepel_0.9.2          fansi_1.0.3            mvtnorm_1.1-3         
##  [16] lubridate_1.9.0        xml2_1.3.3             splines_4.2.2         
##  [19] cachem_1.0.6           knitr_1.41             polyclip_1.10-4       
##  [22] Formula_1.2-4          jsonlite_1.8.4         nloptr_2.0.3          
##  [25] broom_1.0.2            dbplyr_2.2.1           ggforce_0.4.1         
##  [28] compiler_4.2.2         httr_1.4.4             emmeans_1.8.3         
##  [31] backports_1.4.1        assertthat_0.2.1       Matrix_1.5-1          
##  [34] fastmap_1.1.0          gargle_1.2.1           cli_3.6.0             
##  [37] tweenr_2.0.2           htmltools_0.5.4        tools_4.2.2           
##  [40] igraph_1.3.5           coda_0.19-4            gtable_0.3.1          
##  [43] glue_1.6.2             tables_0.9.10          Rcpp_1.0.9            
##  [46] cellranger_1.1.0       jquerylib_0.1.4        vctrs_0.5.1           
##  [49] svglite_2.1.1          nlme_3.1-160           xfun_0.36             
##  [52] lme4_1.1-31            rvest_1.0.3            timechange_0.1.1      
##  [55] lifecycle_1.0.3        googlesheets4_1.0.1    MASS_7.3-58.1         
##  [58] scales_1.2.1           tidygraph_1.2.2        BayesFactor_0.9.12-4.4
##  [61] hms_1.1.2              parallel_4.2.2         TMB_1.9.1             
##  [64] yaml_2.3.6             gridExtra_2.3          pbapply_1.7-0         
##  [67] sass_0.4.4             stringi_1.7.12         highr_0.10            
##  [70] checkmate_2.1.0        boot_1.3-28            rlang_1.0.6           
##  [73] pkgconfig_2.0.3        systemfonts_1.0.4      distributional_0.3.1  
##  [76] evaluate_0.19          lattice_0.20-45        labeling_0.4.2        
##  [79] tidyselect_1.2.0       table1_1.4.3           magrittr_2.0.3        
##  [82] R6_2.5.1               generics_0.1.3         DBI_1.1.3             
##  [85] mgcv_1.8-41            pillar_1.8.1           haven_2.5.1           
##  [88] withr_2.5.0            modelr_0.1.10          crayon_1.5.2          
##  [91] utf8_1.2.2             tzdb_0.3.0             rmarkdown_2.19        
##  [94] viridis_0.6.2          grid_4.2.2             readxl_1.4.1          
##  [97] reprex_2.0.2           digest_0.6.31          webshot_0.5.4         
## [100] numDeriv_2016.8-1.1    munsell_0.5.0          viridisLite_0.4.1     
## [103] bslib_0.4.2